function SAT_SST_Step_04_scale_data(P)

    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_anm_',P.date,'.mat'];
    load(file_load,'raw_SAT_anm')

    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_SST_',P.date,'.mat'];
    load(file_load,'lon','lat','land_fraction')
    
    PP = P;        
    if P.use_qcf_para == 1,  PP.app = 'qcf';  else,  PP.app = 'qcu';  end
    file_load = [SAT_SST_func_IO('data',PP),'/BB98d_parameters_starting_year_',num2str(PP.yr_learn(1)),'_',PP.date,'.mat'];
    load(file_load,'Parameters');

    % *********************************************************************
    % Loop over stations and rescale the data
    % *********************************************************************
    Scl_AT = nan(size(raw_SAT_anm));
    Scl_AT_cowtan = nan(size(raw_SAT_anm));
    
    for ct     = 1:size(raw_SAT_anm,1)
        
        if rem(ct,200) == 0, disp(num2str(ct)); end
        
        % Evaluate the BB98d model ----------------------------------------
        AT = squeeze(raw_SAT_anm(ct,:,:));
        
        if nnz(~isnan(AT)) > 120

            % Interpolate data to run the fitted BB98d model ------------------
            Out_station = pre_process_data_scale(AT);
            
            % Prepare for high-resolution forcing files -----------------------
            clear('T_a','dT_a_dt')
            T_a     = BB98d_forcing_month2high_reso(Out_station.AT_itp);
            T_a     = T_a - nanmean(T_a(1:100));
            dT_a_dt = [(T_a(2:end) - T_a(1:end-1)) ./ (86400/2); 0];
            
            % Evaluate the BB98d model ----------------------------------------
            ct_x = discretize(lon(ct),0:10:360);
            ct_y = discretize(lat(ct),-90:10:90);
            x_fit = squeeze(Parameters(ct_x,ct_y,:));
            P_fit = BB98d_predict(T_a,dT_a_dt,x_fit);
            % TODO: to be extended to grid without a fit
            
            % Mask data to be consistent with the coverage of obs -------------
            P_fit.T_o_m(isnan(Out_station.AT)) = nan;
            
            % Put things into a common file according to years ----------------
            temp = reshape(P_fit.T_o_m,12,numel(Out_station.yrs));
            Scl_AT(ct,:,Out_station.yrs-1849) = temp;
            
            % Do the simple Cowtan scaling ------------------------------------
            Scl_AT_cowtan(ct,:,:) = AT .* (0.86 - 0.25 * land_fraction(ct) / 100);
        end
    end

    file_save = [SAT_SST_func_IO('data',P),'/Scaled_SAT_full_',P.date,'.mat'];
    save(file_save,'Scl_AT','Scl_AT_cowtan','-v7.3');
end

function Out_station = pre_process_data_scale(AT)

    AT  = AT(:);

    % Subset years where both AT and SST has values ---------------
    temp = reshape(AT,12,numel(AT)/12);
    l_yr = any(~isnan(temp),1);
    l_yr(find(l_yr,1):find(l_yr,1,'last')) = 1;
    yrs  = 1850:2020;
    Out_station.yrs  = yrs(l_yr);

    % Linear interpolate if missing values ------------------------
    l_yr = repmat(l_yr,12,1);
    Out_station.AT   = AT(l_yr);

    Out_station.AT_itp  = CDC_interp1(Out_station.AT,[],1);

end
